{ Coccolith Biometrics macros - a set of  macros for NIH-Image developed by Jeremy Young (Natural History Museum, London),  with much assistance from Michal Kucera (Geology Dept., Prague, Goteborg Univs.) and Hsiao-Wen Chung (Bioegineering, Pennsylvania Univ.). Purpose measuring Emiliania huxleyi coccoliths from cross-polarised light images. Queries to jy@nhm.ic.ac.uk. Some macros require Image version 1.54.
* After macro name indicates that operation is influenced by cursor position. 
** Indicates operation will be repeated on every tile of mosaic if cursor is in top left corner (red square). 
Revision 1.0b19, May 1994
}


var {Global Variables}
z:integer {main array reference} 
x,y,x1,x2,x3,x4,y1,y2,y3,y4:real;            {x-y pixel locations, used in various contexts}
cx,cy,xpos,ypos:real;               {x-y pixel locations, used for ellipse center and reporting proc results}
tcount,fcount,pcount:integer;              {counters used by auto-measure macro}
lw,pitdist,oldpitdist,pitval:real;    {variables for proclpit and procpit} 
count,ppv:integer
min,max:real;                                {variables for plotdata}
N,n1,n2,c:integer;                               {counters used in various contexts}
xMos,yMos,Nmos,xN,yN:integer;      {mosaic variables - tile x/y centre, ref no., x,y ref nos.} 
qn1,qn2:integer                                   {tile ref from procsequence};
{N.B. If it is not a mosaic image then xMos/yMos is a reference point set by the cursor location}
greyval:integer;               {used to record greyval of overwritten pixel}
edgethreshold:real;          {threshold ratio of min:(max-min) used by edge finding routine} 
radius:real;                     {radius used by proclpit, must be > max likely c.area length - microns} 
pt:integer                       {ca rim depth threshold - greyscale value}
tl:real;                            {length to test for procedge - microns}
maxrimcaratio:real;        {ratio of local ca radius to searched length for procedge}
cut:real;                          {threshold deviation from mean for rejecting edge points}
hscale,vscale,par:real;          {x and y spatial scales, pixel aspect ratio}
edgeL: real;                     {results of procedge - length (µm)}
pi: real;                           {constant pi, set by procinit}
xc,yc:real;                      {coordinate of ellipse center - microns}
xl,yl,xi,yi:real                 {image & local co-ords - pixels}
x8,y8,v:real                    {results from martin edge}
axisa,axisb:real;             {semi-axes of the ellipse - microns}
theta,k:real;                    {inclination of long axis - radians, clockwise from horiz.; k-various uses}
sin2th,sinth,costh:real;   {temp storage for trig fn values}
found,state,init,mosaic,restCA:boolean;

Procedure ShowParams
Begin
ShowMessage(1:0:0,pt:6:0,' rim depth threshold (pt - 15)','\',2:0:0,radius:6:2,' max c.area l. (radius - 3.0µm)','\',3:0:0,edgethreshold:6:2,' edge threshold (0.1)','\',4:0:0,maxrimcaratio:6:2,' max rim:ca ratio (4)','\',5:0:0,cut:6:2,' edge point elim. (cut - 1.15)','\',0:0:0,' CANCEL','\','\','Defaults are set by procinit - change in macro file')};
End

Procedure procinit;
{function - sets a series of default parameter values} 
{mag, hscale & vscale defaults must be adjusted for your microscope}
Var
mag:integer;
Begin
  par:=1.05 {Pixel Aspect Ratio on my system};
  mag:=GetNumber('What is magnif (objxopt)?',160);
  hscale:=mag*0.0964;  {usual value on my microscope is 15.424pixels/micron}
  vscale:=hscale/par;   {usual value on my microscope is 14.689} 
  init:=1;
  z:=0;       {actual value of z used will usually be set later}
  SetCounter(30);
  SetFontSize(9);
  SetText('Left Justified With Background');
  SetBackGroundColor(255);
  SetOptions('X-Y Center Length Angle Min/Max User1');
  SetUser1Label('RimW');
  SetLineWidth(1);
  pi:=3.141592654;
  pt:=15;
  maxrimcaratio:=4;  
  cut:=1.15;                    
  radius:=3.0;
  edgethreshold:=0.05;
  ShowParams;
end;


{====================Mosaic Macros=============================}

procedure procn2xy(Nmos);
{given the ref number Nmos finds the x,y loc of the centre of the mosaic tile}
Begin
yN:=round((Nmos-3.5)/6);
xN:=Nmos-(yN*6)-1;
xMos:=55+xN*105;
yMos:=55+yN*105;
end;

procedure procxy2n(x,y);
{given an x,y loc finds the ref number Nmos of the mosaic tile}
Begin
xN:=round((x-55)/105);
yN:=round((y-55)/105);
Nmos:=xN+6*yN+1;
xMos:=55 +xN*105;
yMos:=55+yN*105;
end;

Procedure procSetSequence(x,y);
{function  - determine which specimen(s) to work on}
var w,h:integer;
Begin
GetPicSize(w,h);
IF (w*1000+h)=636530 THEN BEGIN {test if this a mosaic image}
    mosaic:=1; 
    IF (x<20) AND (y<20) THEN BEGIN  {If cursor in top left corner work on all specs}
       qn1:=1; qn2:=30;
       xMos:=55; yMos:=55;
    END ELSE BEGIN                            {If cursor is elsewhere  work on spec in tile}
       procxy2n(x,y);
       qn1:=nMos; qn2:=nMos;
    END;
END ELSE BEGIN                                {If not a mosaic work on spec under cursor} 
   mosaic:=0;
   qn1:=z+1; qn2:=z+1;
   xMos:=x; yMos:=y;
END;
End;

Macro 'TAKE MOSAIC PHOTO * [3]' 
{function - take "photo" of specimen and save to mosaic} 
Begin
  StopCapturing;
  IF init=0 then procinit;  
  KillRoi;
  GetMouse(x,y);
  n:=nMos;
  ProcSetSequence(x,y);
  nMos:=n;  
  MakeRoi(xMos-50,yMos-50,100,100);
  Copy;
  NextWindow;
  procn2xy(nMos); 
  MakeRoi(xMos-50,yMos-50,100,100);
  Paste;
  MoveTo(xMos-50,yMos+50);
  setForegroundColor(2);  
 Write(Nmos:0:0);  
  Nmos:=Nmos+1;
   If Nmos>30 then begin beep; SaveAs; End;
  NextWindow;
  ShowMessage(Nmos,'-Nmos');
  KillRoi;
  {StartCapturing;}
End

macro 'Make blank mosaic window [4]'
{function - set up blank mosaic window}
Begin
   StopCapturing;
   SetNewSize(635,530);
   SetBackGroundColor(255);
   MakeNewWindow('Mosaic');
   SetForeGroundColor(1);
   MakeRoi(0,0,10,10);Fill; {Make red corner, if mouse is here sequencer will run}  
   KillRoi;   
   Nmos:=1;
   procn2xy(Nmos); 
   PutMessage('If window is not scaled then SetMag {f3}. Close Other windows before taking mosaic photos');
   SaveAs;
   NextWindow;
End;

Macro 'Set number of next tile * [N]'
{function - set nMos, useful e.g. when compiling mosaics over a period of time} 
Begin
GetMouse(x,y);
procxy2n(x,y); 
Nmos:=GetNumber('tile ref (1-30)',Nmos);
procn2xy(Nmos);
greyval:=GetPixel(xMos,yMos);
PutPixel(xMos,yMos,1);
ShowMessage(Nmos,'=Nmos','\',yN,xN,'   yN,xN','\', xMos, yMos);
WaitForTrigger;
PutPixel(xMos,yMos,greyval);
end;

macro 'Retrieve tile from duplicate * [D]' 
{func
tion - retrieve a duplicate tile from another copy of the mosaic} 
Begin
  GetMouse(x,y);
  procxy2n(x,y);
  NextWindow;
  MakeRoi(xMos-50,yMos-50,100,100);
 Copy;
  NextWindow;
  procn2xy(Nmos); 
  MakeRoi(xMos-50,yMos-50,100,100);
  Paste;
  KillRoi;
End

macro 'Centre mosaic specimen * [C]'
{function - correct position  of off-centre specimens}
Begin
GetMouse(x,y);
procxy2n(x,y);
MakeRoi(x-40,y-40,80,80);
Copy;
MakeRoi(xMos-40,yMos-40,80,80);
Paste;
  MoveTo(xMos-50,yMos+50);
  SetForegroundColor(2);  
 Write(Nmos:0:0);  
End;

Macro 'Show tile centres** [Q]'
Begin
GetMouse(x,y);
ProcSetSequence(x,y);
FOR n:=qn1 TO qn2 DO BEGIN
  PutPixel(xMos,yMos,2);
  procn2xy(n+1);
END;
END;

Macro'(-'; 

{======== ========== Biometric Macros ============= ============}


Procedure procEllipseDraw(axisa,axisb,cx,cy,alpha,n,n2,colour:real,restore:boolean);
var
ang:real;
Begin
FOR count:=1 to N do begin
 ang:=2*pi*(count/n);
 xl:=axisa*cos(ang);
 yl:=axisb*sin(ang);
 x:=cx+cos(ALPHA)*xl - sin(ALPHA)*yl;
 y:=cy+sin(ALPHA)*xl + cos(ALPHA)*yl;
 IF NOT Restore THEN BEGIN
   greyval:=rUser1[n2+count];
   PutPixel(x,y,greyval);
 END ELSE BEGIN
   rUser1[n2+count]:=GetPixel(x,y);
   PutPixel(x,y,colour);
 End;
end;
End;

procedure procMartinEdge(a,b,x0,y0:real); 
var aSQ,bSQ,e,k,kSQ,q,xs,ys,k2,q2,k3,q3,q4:real;
{x8,y8, position of nearest point on ellipse (a,b) to x0,y0. v distance between x0,y0 & x8,y8} 
begin
 aSQ:=a*a; bSQ:=b*b;
 x1:= a*b*x0/sqrt(bSQ*x0*x0+aSQ*y0*y0);
 y1:= y0*x1/x0;
 e:= sqrt(aSQ-bSQ);
 k:= y0/(x0-e*e/a); kSQ:=k*k;
 q:= -k*e*e/a;
 x2:= (-k*q*aSQ-abs(q)/q*a*b*sqrt(bSQ-q*q+aSQ*kSQ))/(bSQ+aSQ*kSQ);
 y2 := k*x2+q;
 k:= (y2-y1)/(x2-x1); kSQ:=k*k;
 q:= -k*((y1+y2)/2)-(x1+x2)/2;
 y3:= (-k*q*bSQ+a*b*sqrt(aSQ-q*q+bSQ*kSQ))/(aSQ+bSQ*kSQ);
 x3:= -k*y3-q;
 k2:= (y2-y3)/(x2-x3);
 k3:= (y1-y3)/(x1-x3);
 q2:= -k2*(y2+y3)/2-(x2+x3)/2;
 q3:= -k3*(y1+y3)/2-(x1+x3)/2;
 ys:= (q2-q3)/(k3-k2);
 xs:= -k2*ys-q2;
 k:= (ys-y0)/(xs-x0); kSQ:=k*k;
 q4:= y0-k*x0;
 x8:= (-k*q4*aSQ+a*b*sqrt(bSQ-q4*q4+aSQ*kSQ))/(bSQ+aSQ*kSQ);
 y8:= k*x8+q4;
 v:= sqrt(sqr(x8-x0)+sqr(y8-y0));
end;

procedure procImageToLocal(xq,yq,alpha,cx,cy);
Begin
xl:=cos(ALPHA)*(xq-cx) - sin(ALPHA)*(yq-cy);
yl:=sin(ALPHA)*(xq-cx) + cos(ALPHA)*(yq-cy);
End;

procedure procLocalToImage(xl,yl,alpha,cx,cy);
Begin
xi:=cos(ALPHA)*xl + sin(ALPHA)*yl + cx;
yi:=-sin(ALPHA)*xl + cos(ALPHA)*yl + cy;
End;


Procedure procDrawAxes(z1);
Begin
   SetForegroundColor(1);  
   MoveTo(rX[131],rY[131]);
   LineTo(rX[132],rY[132]);               {display ca length, as red line}
   SetForegroundColor(2); 
   MoveTo(rX[133],rY[133]);
   LineTo(rX[134],rY[134]);               {display ca width, as green line}
  x1:=rX[131]; y1:=rY[131];
  SetForeGroundColor(5);
  theta:=rAngle[z1]*pi/180;
  x2:=x1-rUser1[z1]*cos(theta)*hscale;
  y2:=y1-rUser1[z1]*sin(theta)*vscale;
  MoveTo(x1,y1); LineTo(x2,y2);     {display rimwidth, as pale blue line}
End;


Procedure procTogglepixelvals;
{for use after running main macro, to show/hide values}
Begin
 IF state=1 then  begin
  FOR n:=tcount+1 to pcount do begin
    PutPixel(rX[n],rY[n],rLength[n]);
  end; 
  FOR n:=101 to tcount do begin
    PutPixel(rX[n],rY[n],rLength[n]);
   state:=0;
  end;
 end else begin
  FOR n:=tcount+1 to pcount do begin
    PutPixel(rX[n],rY[n],2);
    end; 
    FOR n:=101 to tcount do begin
       PutPixel(rX[n],rY[n],1);
    end;
    state:=1
 end; 
END;

procedure procedge(xt1,yt1,xt2,yt2);
{function - find edge of coccolith shield}
{Calculates a threshold value based on the max & min values and looks for series of pixels above it. N.B. Works from outside inwards, i.e. from background toward specimen}
var
trip,thresh,xtn1,ytn1:integer
Begin
RequiresVersion(1.54);
IF mosaic THEN BEGIN   {truncate line if it extends out of tile}
IF xt1-xMos>50 THEN BEGIN xtn1:=xMos+49; yt1:=(xtn1-xt2)/(xt1-xt2)*(yt1-yt2)+yt2;xt1:=xtn1;END;
IF xt1-xMos<-50 THEN BEGIN xtn1:=xMos-49; yt1:=(xtn1-xt2)/(xt1-xt2)*(yt1-yt2)+yt2;xt1:=xtn1;END;
IF yt1-yMos>50 THEN BEGIN ytn1:=yMos+49; xt1:=(ytn1-yt2)/(yt1-yt2)*(xt1-xt2)+xt2;yt1:=ytn1;END;
IF yt1-yMos<-50 THEN BEGIN ytn1:=yMos-49; xt1:=(ytn1-yt2)/(yt1-yt2)*(xt1-xt2)+xt2;yt1:=ytn1;END;
END;
PutPixel(xt1,yt1,3);
MakeLineRoi(xt1,yt1,xt2,yt2);
GetPlotData(count,ppv,min,max);
c:=-1;
found:=0;
REPEAT
   c:=c+1;
UNTIL PlotData[c]=max;       {advance to max value}
If c=0 THEN c:=1;
thresh:= max-(max-min)*edgethreshold;
REPEAT
   c:=c+1;
   IF PlotData[c]<(thresh-2) THEN BEGIN
       trip:=trip+1;
   END ELSE BEGIN
       trip:=0;
   END;
UNTIL (trip=4) OR (c>=count);
IF trip=4 THEN BEGIN
   c:=c-4;
   found:=1;
   xpos:=xt1+(c/count*(xt2-xt1));
   ypos:=yt1+(c/count*(yt2-yt1));
   putpixel(xpos,ypos,1);
   procImagetoLocal(xpos,ypos,-theta,cx,cy);
   xl:=abs(xl);  yl:=abs(yl);
   procMartinEdge(axisa*hscale,axisb*hscale,xl,yl);
   edgeL:=v;
end;
end

procedure procpit(x1,y1,n:integer);
{finds lowest val in (2n+1)x(2n+1) box, jumps to it, repeats until at local inverse peak, colours it red}
  var
pkx,pky,oldval,xi,yi,val: integer;
Begin
  pitval:=255;
REPEAT
 oldval:=pitval; 
 for xi:=x1-n to x1+n DO begin
   for yi:=y1-n to y1+n DO begin
     val:=GetPixel(xi,yi);
     IF val<pitval THEN begin
       pitval:=val;
       pkx:=xi;
       pky:=yi;
     end;
   end;
 end;        {end of nest of loops finding lowest value in box}
 x1:=pkx;
 y1:=pky;
UNTIL pitval=oldval;
xpos:=x1;
ypos:=y1;
END 

Procedure ProcLMinimum (x1,y1,x2,y2); 
{function - find xy loc of minimum value on line}
var
Begin
RequiresVersion(1.54);
MakeLineRoi(x1,y1,x2,y2);
GetPlotData(count,ppv,min,max);
c:=1;
REPEAT
   c:=c+1;
UNTIL  PlotData[c]=min;
ypos:=y1+(c/count*(y2-y1));
xpos:=x1+(c/count*(x2-x1));
end;

Procedure proclpit(x1,y1,x2,y2,pt1); 
{function - find 1st major pit on line}
var
pitc,peaksep,nn,oldval,val:integer;
Begin
RequiresVersion(1.54);
MakeLineRoi(x1,y1,x2,y2);
GetPlotData(count,ppv,min,max);
  pitval:=PlotData[0];
  peaksep:=15;
  n2:=4;
  pitc:=0;
  {pitc is the lowest value found - c is the value currently being checked}; 
  nn:=0;
REPEAT
   c:=pitc+nn; 
   val:=PlotData[c];
   nn:=nn+1;
   IF val>(pitval+pt1) THEN peaksep:=n2  else IF val<pitval THEN begin
      pitval:=val;
      pitc:=c; 
      nn:=1; 
   end;
UNTIL (c=count) OR (nn>peaksep);
IF (c<count) AND (peaksep=n2) THEN Begin
 ypos:=y1+(pitc/count*(y2-y1));
 xpos:=x1+(pitc/count*(x2-x1));
 pitdist:=SQRT(sqr(x1-xpos)+sqr(y1-ypos)); 
 found:=1;
End;
End;

procedure procEllipseFitHW;
{function - fit ellipse to a ring of points, given the location of the centre}
{N.B. Points come from proclpit & procpit, centre from analysis of pits} 
{pixel co-ordinates of the centre and boundary points are stored in arrays rX and rY, and are converted into micron co-ordinates with the centre as the origin in arrays rUser1 & rUser2}

var
  i,j:integer;                 {array index general}
  a11,a12,a13:real;      {for matrix A}
  a21,a22,a23:real;
  a31,a32,a33:real;
  det:real;                     {determinant of matrix A}
  b1,b2,b3:real;            {for vector B}
  i11,i12,i13:real;        {for matrix A^(-1)}
  i21,i22,i23:real;
  i31,i32,i33:real;
  c1,c2,c3:real;           {for coefficients}

begin                 {set coordinate system}
  InvertY(false);
  xc:=rX[z]/hscale;
  yc:=rY[z]/vscale;      {coordinates for center point in microns}
  for i:=101 to pCount do begin
    rUser1[i]:=rX[i]/hscale-xc;
    rUser2[i]:=rY[i]/vscale-yc;         {coordinate translation, pixels -> microns wrt ellipse centre}
  end;
  
  for i:=101 to pCount do begin
    a11:=a11+rUser1[i]*rUser1[i]*rUser1[i]*rUser1[i];
    a12:=a12+rUser1[i]*rUser1[i]*rUser1[i]*rUser2[i];
    a13:=a13+rUser1[i]*rUser1[i]*rUser2[i]*rUser2[i];
    a23:=a23+rUser1[i]*rUser2[i]*rUser2[i]*rUser2[i];
    a33:=a33+rUser2[i]*rUser2[i]*rUser2[i]*rUser2[i];
    b1:=b1+rUser1[i]*rUser1[i];
    b2:=b2+rUser1[i]*rUser2[i];
    b3:=b3+rUser2[i]*rUser2[i];        {calculation of matrices A & B}
  end;
  a21:=a12;
  a22:=a13;
  a31:=a13;
  a32:=a23;

  {determinant of A}
  det:=a11*a22*a33+a21*a32*a13+a31*a12*a23
      -a11*a32*a23-a21*a33*a12-a31*a22*a13;
  if det=0 then begin
    PutMessage('Sorry.  Singular matrix found');
    Exit;
  end;
  i11:=(a22*a33-a32*a23)/det;
  i12:=(a32*a13-a12*a33)/det;
  i13:=(a12*a23-a22*a13)/det;        {inversion of matrix A}
  i22:=(a11*a33-a31*a13)/det;
  i23:=(a21*a13-a11*a23)/det;
  i33:=(a11*a22-a21*a12)/det;
  i21:=i12;
  i31:=i13;
  i32:=i23;

  c1:=i11*b1+i12*b2+i13*b3;
  c2:=i21*b1+i22*b2+i23*b3;          {matrix multiplication to get coefficients}
  c3:=i31*b1+i32*b2+i33*b3;

  if c2=0 then theta:=0            {find orientation in real world, between 145!}
  else if c1=c3 then theta:=0.7854      {45!}
  else theta:=arctan(c2/(c1-c3))/2;

  {find length and width in microns}
  if theta=0 then begin
    axisa:=1/sqrt(c1);
    axisb:=1/sqrt(c2);
  end
  else begin
    sin2th:=sin(2*theta);
    axisa:=1/sqrt((c1+c3+c2/sin2th)/2);
    axisb:=1/sqrt((c1+c3-c2/sin2th)/2);
  end;
  if axisa<axisb then begin
     k:=2*axisb;
     axisb:=axisa;
     axisa:=k/2;
     theta:=theta+pi/2;
  end;
  rMax[z]:=2*axisa;
  rMin[z]:=2*axisb;
  rAngle[z]:=180*theta/pi;           {rAngle = degrees clockwise from horizontal of the long axis}
end;

Procedure procCalcPixelLocs(z1)
Begin
  {calculate ellipse co-ordinates and store in array(as pixel values)}
  theta:= rAngle[z1]*pi/180;
  costh:=cos(theta);
  sinth:=sin(theta);
  xc:=rX[z1]/hscale;
  yc:=rY[z1]/vscale;                {coordinates for center point in microns}
  axisa:=rMax[z1]/2;
  axisb:=rMin[z1]/2;
  rX[131]:=(xc-axisa*costh)*hscale;
  rY[131]:=(yc-axisa*sinth)*vscale;
  rX[132]:=(xc+axisa*costh)*hscale;
  rY[132]:=(yc+axisa*sinth)*vscale;
  rX[133]:=(xc-axisb*sinth)*hscale;
  rY[133]:=(yc+axisb*costh)*vscale;
  rX[134]:=(xc+axisb*sinth)*hscale;
  rY[134]:=(yc-axisb*costh)*vscale;
end;

Procedure procCalcRimWidth(z,n1,n2);
 {determines a sensible rim width, from a set of local estimates}
Var
na,na2,ave,ave2:real;
Begin   
FOR n:=n1 TO n2 DO BEGIN
     IF rlength[n]<rMax[z] then begin
       ave:=ave+rlength[n];
       na:=na+1;
     END;
   END;
   ave:=ave/na;
   FOR n:=n1 TO n2 DO BEGIN
     IF (rlength[n]>ave/cut) AND (rlength[n]<ave*cut) then begin
        PutPixel(rMin[n],rMax[n],4);      {colour pixels within cut, yellow} 
        ave2:=ave2+rlength[n];
        na2:=na2+1;
      END ELSE PutPixel(rMin[n],rMax[n],1);    {colour pixels outside cut, red} 
   END;
   IF na2<2 THEN BEGIN
     na2:=0; ave2:=0;
     FOR n:=n1 TO n2 DO BEGIN
       IF rLength[n]>ave THEN BEGIN
          PutPixel(rMin[n],rMax[n],2);    {colour high edgeL pixels - as used - green} 
          ave2:=ave2+rlength[n];
          na2:=na2+1;
        END ELSE PutPixel(rMin[n],rMax[n],1);
     END;
   END; 
   rUser1[z]:=ave2/na2;
   rLength[z]:=rMax[z]+2*rUser1[z];
 END;

Procedure procLineError;
Begin
PutMessage('Select line of correct length before running macro');
Exit; END;
IF x1=-1 THEN procLineError;


Macro 'AUTO-MEASURE** [5]'
{function - automatically measure central area length and width of coccolith, rim width, and orientation. Imput a single pint inside the central area of the cocolith}
{rX,rY,rLength,tcount -x,y,pixel vals, counter -  of the pits
 rUser1,rUser2,rMax,fcount -  x,y,pixel vals, counter - of points on the ca edge
 before doing ellipse fit the ca edge points are appended to the list of pit points and the counter pcount applies to entire series.
N.B. tcount, fcount & pcount all run from 101 up since they form indices to a work area of the measurement arrays.}

var
ang,Dang:real;                                           {angle and Æangle for linepit procedure} 
loc,loc2,nloc,pixelradius,diag,l1,l2,lax,lax1,lax2:integer;
BEGIN
RequiresVersion(1.54);
IF init=0 then procinit;
IF rCount<149 THEN BEGIN SetCounter(150); z:=0;end;
pixelradius:=radius*hscale;
Getmouse(x1,y1);
procSetSequence(x1,y1);
FOR z:=qn1 TO qn2 DO BEGIN                           {this causes repeat of full sequence if working on mosaic}
oldpitdist:=pixelradius;
ang:=0;
rX[101]:=xMos;
rY[101]:=yMos;
rLength[101]:=GetPixel(xMos,yMos);
loc:=xMos*1000+yMos;
loc2:=0;
tcount:=101; 
fcount:=101;
Dang:=(20/360)*(2*pi);
{pixelradius is the radius of the area/line length searched for pit value};
REPEAT
  found:=0;
  proclpit(xMos,yMos,xMos+cos(ang)*pixelradius,yMos+sin(ang)*pixelradius, pt);
  IF (found=1) AND (pitdist<2*oldpitdist) THEN begin
     rAngle[fcount]:=360/(2*pi)*ang;
     ang:=ang + Dang;            {this causes the next angle to be skipped}
     oldpitdist:=pitdist;
     fcount:=fcount+1;
     rUser1[fcount]:=xpos;
     rUser2[fcount]:=ypos;
     rMax[fcount]:=pitval;
     procpit(xpos,ypos,2);
     nloc:=xpos*1000+ypos;
    IF((nloc<>loc)  AND (nloc<>loc2))THEN begin  {checks that this is a new pit}
       tcount:=tcount+1;
       loc:=nloc;
       rX[tcount]:=xpos;
       rY[tcount]:=ypos;
       rLength[tcount]:=pitval;
    END; 
    loc2:=rX[102]*1000+rY[102];
  END;
  ang:=ang+Dang;
UNTIL  ang>=2*pi;
{
Next routine, finds the two most widely separated peaks and centre of line between them
}
IF tcount>103 THEN BEGIN
l2:=0;
For n1:=102 to tcount do begin
  For n2:=(n1+1) to tcount do begin
    l1:=sqrt(sqr(rX[n1]-rX[n2])+sqr(rY[n1]-rY[n2]));
    IF l1>l2 then begin
      l2:=l1;
      lax1:=n1;
      lax2:=n2;
    End;
  End;
End;
cx:=0.5*(rX[lax1]+rX[lax2]);
cy:=0.5*(rY[lax1]+rY[lax2]);
rX[z]:=cx;
rY[z]:=cy;
State:=0;
For n:= 102 to fcount do begin
   rX[tcount-101+n]:=rUser1[n];
   rY[tcount-101+n]:=rUser2[n];
   rLength[tcount-101+n]:=rMax[n];
end;
pcount:=tcount-101+n;
killroi;
IF (pcount>106)  THEN BEGIN       {With <5 points procEllipseFit won't work}
    procEllipseFitHW;

{next set of routines finds rim width}
    procCalcPixelLocs(z);
   For n:=102 to tcount DO BEGIN
      x1:=rX[n]; y1:=ry[n];
      x2:=(x1-cx)*maxrimcaratio+cx; 
      y2:=(y1-cy)*maxrimcaratio+cy;
      procedge(x2,y2,x1,y1);
      rLength[n+29]:=edgeL/hscale;
      rMin[n+29]:=xpos; rMax[n+29]:=ypos;
    End;
  procCalcRimWidth(z,131,tcount+29);
 
{next set of routines displays results}
   SetForegroundColor(255);  
   MoveTo(rX[lax1],rY[lax1]);
   LineTo(rX[lax2],rY[lax2]);                 {display diagonal betwen pits, as black line}
   procDrawAxes(z);
   proctogglepixelvals;
 END;                          {end of routines dependant on finding ellipse axes}
 END;                          {end of routines dependant on finding ellipse centre}
procn2xy(z+1);
END;                           {end of mosaic-loop}
IF n2=30 then beep;
END;
end;


Macro'SEMI-MANUAL MEASURE * [6]'
{N.B. Input is a line region of interest - a line going through the ellipse axis- and a point inside the coccolith ring indicated by the centre of the line selection - i.e. centre of line selection must be inside the c.area. Output results contain central area parameters and the coccolith rim width}

var
  dx,dy,x0,y0:real;
begin
IF init=0 then procinit;  
found:=0;
getline(x1,y1,x2,y2,lw);
IF x1=-1 THEN procLineError;
procxy2n(x1,y1);
z:=nMos;                       {selects correct spec if working on mosaic}
rAngle[z]:=arctan((y1-y2)/(x1-x2))*180/pi;
x0:=0.5*(x1+x2);
y0:=0.5*(y1+y2);
proclminimum(x0,y0,x1,y1);
rX[131]:=xpos; rY[131]:=ypos;
proclminimum(x0,y0,x2,y2);
rX[132]:=xpos; rY[132]:=ypos;

{this has found the position of two maxima on the original line
  Next routine finds mid point between maxima, and new maxima on perpendicular line}

rMax[z]:=sqrt(sqr(rX[131]-rX[132])+sqr(PAR*(rY[131]-rY[132])))/15.424;
cx:=(rX[131]+rX[132])/2;
cy:=(rY[131]+rY[132])/2;
rX[z]:=cx; rY[z]:=cy;

dx:=(x1-x2)*.5;  
dy:=(y1-y2)*.5;
x3:=cx-dy;   y3:=cy+dx;
x4:=cx+dy;  y4:=cy-dx;
proclminimum(cx,cy,x3,y3);
rX[133]:=xpos; rY[133]:=ypos;
proclminimum(cx,cy,x4,y4);
rX[134]:=xpos; rY[134]:=ypos;
rMin[z]:=sqrt(sqr(rX[133]-rX[134])+sqr(PAR*(rY[133]-rY[134])))/hscale;
rUser1[z]:=0;
ProcCalcPixelVals(z);
procDrawAxes(z);
end;


Macro 'Correct long axis * [7]';
{function - manually input/correct measurements associated with the long axis} 
Begin
IF init=0 then procinit;  
GetLine(x1,y1,x2,y2,lw);
IF x1=-1 THEN procLineError;
procxy2n(x1,y1);
rMax[nMos]:=sqrt(sqr(x1-x2)+sqr(PAR*(y1-y2)))/hscale;
rX[nMos]:=0.5*(x1+x2);
rY[nMos]:=0.5*(y1+y2);
rAngle[nMos]:=arctan((y1-y2)/(x1-x2))*180/pi;
SetForeGroundColor(1); Fill;
KillRoi;
End;

Macro 'Correct short axis * [8]';
{function - manually input/correct short axis length}
Begin
IF init=0 then procinit;  
GetLine(x1,y1,x2,y2,lw);
IF x1=-1 THEN procLineError;
procxy2n(x1,y1);
rMin[nMos]:=sqrt(sqr(x1-x2)+sqr(PAR*(y1-y2)))/hscale;
SetForeGroundColor(2); Fill;
KillRoi;
End;

Macro 'Correct Rim Width * [9]'
{function - manually input rim width, needs central area measurements to have been made already}
Begin
GetMouse(x,y);
Procxy2n(x,y);
z:=nMos;
axisa:=0.5*rMax[z]*hscale; axisb:=0.5*rMin[z]*hscale;
theta:=rAngle[z]/360*2*pi;
cx:=rX[z]; cy:=rY[z];
v:=rUser1[z]*hscale;
  putpixel(x,y,5);
  procImagetoLocal(x,y,-theta,cx,cy);
  x:=abs(xl); y:=abs(yl);
  procMartinEdge(axisa,axisb,x,y); 
  procEllipseDraw(axisa+v,axisb+v,cx,cy,theta,30,250,4,1);
  rUser1[z]:=v/hscale;
  rLength[z]:=rMax[z]+2*rUser1[z];
End;


Macro 'Draw Axes** [A]'
Begin
IF NOT init THEN procinit;
GetMouse(x,y);
ProcSetSequence(x,y);
FOR z:=qn1 TO qn2 DO BEGIN
  ProcCalcPixelVals(z);
  ProcDrawAxes(z);
  procn2xy(z+1);
END;
END;


Macro 'DrawCaEllipse * [E]'
Begin
GetMouse(x,y);
Procxy2n(x,y);
z:=nMos;
axisa:=0.5*rMax[z]*hscale; axisb:=0.5*rMin[z]*hscale;
theta:=rAngle[z]/360*2*pi;
cx:=rX[z]; cy:=rY[z];
restCA:=NOT restCA;
procEllipseDraw(axisa,axisb,cx,cy,theta,30,200,5,restCA);
End;


Macro 'Toggle pixel vals [T]'
{displays/hides the pixels found by proclpit and procpit on the last measured specimen}
Begin
RequiresVersion(1.54);
procTogglePixelVals;
END;

Macro 'Surface plot specimen * [S]'
{function - produce surface plot of a single specimen}
Begin
RequiresVersion(1.50);
GetMouse(x,y);
n:=picnumber;
GetRoi(x1,y1,x2,y2);
IF x2<10 THEN BEGIN 
  ProcSetSequence(x,y);
  MakeRoi(xMos-48,yMos-48,96,90);
END;
SetNewSize(400,400);
Invert;
SurfacePlot;
ShowMessage('Press cmnd-. while clicking to keep surface plot');
WaitForTrigger;
n2:=picnumber;
ChoosePic(n);
Invert;
KillRoi;
ChoosePic(n2);
Dispose;
ChoosePic(n);
Wait(0.2)
End;


 
Macro 'Set parameters [P]'
{function - to allow operator to change key parameters/tune operation for different species or samples}
var 
i: integer;
Begin
ShowParams;
i:=GetNumber('Choose Param - values window',0);
IF i=1 THEN BEGIN 
ShowMessage('pt is depth (in grey values) of the pit/trough that proclpit looks for','\','Set LOW (<10) to analyze small specimens','\','Set HIGH (>20) to analyze messy or bright specimens','\',15,'=default','\',pt,'=current value');
pt:=GetNumber('rim depth threshold (pt - eg 15)',pt);
END;
IF i=2 THEN BEGIN 
ShowMessage('radius is the distance (µm) from reference point searched by proclpit. It should  Åmax likely c.area length');
radius:=GetNumber('max c.area l. (radius - eg 3.0µm',radius);
END;
IF i=3 THEN BEGIN
ShowMessage('edge threshold sets value used to detect c.area edge, = ratio of height of local background and rim values','\',' Set LOW (<0.1) to move edge detection out','\','Set HIGH (>0.1) to move edge detection in'); 
edgethreshold:=GetNumber('edge threshold (eg 0.1)',edgethreshold);
END; 
IF i=4 THEN BEGIN
ShowMessage('max rim:ca ratio','\','This sets the distance from the coccolith centre searched by procEdge, as a mutiplier of the centre-ca distance');
maxrimcaratio:=GetNumber('max rimcaratio(eg  4)',maxrimcaratio);
END;
IF i=5 THEN BEGIN
ShowMessage('cut is the criterion used to eliminate edge points');
cut:=GetNumber('edge point elim. value (cut - eg 1.15)',cut);
END;
ShowParams;
End;

Macro 'Initialise [I]';
Begin
procinit;
End;

Macro '(-'

{==============Utility macros=====================================}

macro 'Kill ROI [f2]'  begin KillRoi end;

macro 'Magn + ScaleBar [f3]'
{function - set spatial scale, based on microscope lenses used}
   var
obj,width,height:integer;
opt,mag:real;
Begin
  GetPicSize (width,height);
  n:=GetNumber('Draw Scale Bar? (1=Yes, 0=No)',0);
If n<>2 then begin 
 obj:=GetNumber('Objective Lens',100);
  opt:=GetNumber('Optovar Lens',1.6);
  mag:=obj*opt;
  hscale:=mag*0.0964;
  vscale:=hscale/PAR; 
  ShowMessage('1.05=Aspect Ratio');
  SetScale(hscale,'µm');
  SetScale(0,'');     {opens dialog box -  needed as Aspect Ratio can't be set by macros};
 End;  
 If  n>=1 then begin
    SetForegroundColor(0);
    SetBackGroundColor(255);
    MakeRoi (10,height-15,10*hscale,5);
    Fill;
    MoveTo(10*hscale+20,height-12);
    Write('10µm');
  end;
end;

Macro 'Shrink ROI [f5]' begin InsetRoi(10) end;

Macro  'Select ROI * [f6]'
  var
x,y: integer;
Begin
  GetMouse(x,y);
  MakeRoi(x-50,y-50,100,100);
end:

Macro 'Grow ROI [f7]' begin InsetRoi(-10) end;

macro 'Reset Greymap[f9]'  begin ResetGrayMap end;

macro 'Open Grey + 4 ramp LUT[f10]';  begin Open('grey+4thinramps'); end;

macro 'Add green line to LUT [f11]'
begin
RedLUT[2]:=0;GreenLUT[2]:=255;BlueLUT[2]:=0;
UpdateLUT;
end;